# A1_political_knowledge_multiple_media.R
# Created: 2022-8-31 Taka-aki Asano
# Last Modified: 2022-8-31

# package
require("psych")
require("magrittr")
require("nnet")
require("coefplot")


# explanatory variables
voter2019$FEMALE <- ifelse(voter2019$SEX == 2, 1, 0)
voter2019$GENERATION <- NA
voter2019$GENERATION[voter2019$AGE <= 29] <- 1
voter2019$GENERATION[voter2019$AGE >= 30 & voter2019$AGE <= 39] <- 2
voter2019$GENERATION[voter2019$AGE >= 40 & voter2019$AGE <= 49] <- 3
voter2019$GENERATION[voter2019$AGE >= 50 & voter2019$AGE <= 59] <- 4
voter2019$GENERATION[voter2019$AGE >= 60 & voter2019$AGE <= 69] <- 5
voter2019$GENERATION[voter2019$AGE >= 70] <- 6
voter2019$EDUCATION[voter2019$EDUCATION == 5] <- NA
voter2019$INTEREST_RE <- recode(voter2019$INTEREST, `1` = 5, `2` = 4, `3` = 3, `4` = 2, `5` = 1)


# correlations between media
cor.test(voter2019$NEWSPAPER, voter2019$TV)
cor.test(voter2019$NEWSPAPER, voter2019$ONLINE)
cor.test(voter2019$TV, voter2019$ONLINE)


# correlations between media choice and education
cor.test(voter2019$EDUCATION, voter2019$NEWSPAPER)
cor.test(voter2019$EDUCATION, voter2019$TV)
cor.test(voter2019$EDUCATION, voter2019$ONLINE)


# correlations between media choice and generation
cor.test(voter2019$GENERATION, voter2019$NEWSPAPER)
cor.test(voter2019$GENERATION, voter2019$TV)
cor.test(voter2019$GENERATION, voter2019$ONLINE)


# correlations between media choice and political interest
cor.test(voter2019$INTEREST_RE, voter2019$NEWSPAPER)
cor.test(voter2019$INTEREST_RE, voter2019$TV)
cor.test(voter2019$INTEREST_RE, voter2019$ONLINE)


# multinomial logistic regression
voter2019_mlogit <- dplyr::select(
  voter2019, Knowledge_Class, FEMALE, GENERATION, 
  EDUCATION, INTEREST_RE, NEWSPAPER, TV, ONLINE
)
voter2019_mlogit <- na.omit(voter2019_mlogit)

## interaction term between Newspaper and TV
multi_logit1 <- multinom(
  Knowledge_Class ~ FEMALE + GENERATION + 
    EDUCATION + INTEREST_RE + NEWSPAPER + TV + NEWSPAPER * TV, 
  data = voter2019_mlogit, reflevel = "Class1"
)
summary(multi_logit1)
multi_logit1_plot <- broom::tidy(multi_logit1, 
                                 conf.int = TRUE, 
                                 exponentiate = FALSE) %>% 
  dplyr::filter(term != "(Intercept)") %>% 
  mutate(term = recode(term, 
                       `FEMALE` = "Female", 
                       `GENERATION` = "Age",  
                       `EDUCATION` = "Education", 
                       `INTEREST_RE` = "Political Interest", 
                       `NEWSPAPER` = "Newspaper", 
                       `TV` = "TV", 
                       `NEWSPAPER:TV` = "Newspaper * TV")) %>% 
  mutate(term = factor(term, 
                       levels = rev(c("Female", "Age", "Education", 
                                      "Political Interest", "Newspaper", 
                                      "TV", "Newspaper * TV")))) %>% 
  ggplot(aes(x = term, y = estimate, 
             ymin = conf.low, ymax = conf.high)) + 
  geom_point() + geom_hline(yintercept = 0, color = "gray", linetype = "dashed") + 
  geom_pointrange() + coord_flip() + facet_wrap(~y.level, ncol = 3) + 
  labs(x = "Explanatory Variables", y = "Estimation", 
       title = "Multinomial Logistic Regression \nBaseline: Class1") + 
  theme(plot.title = element_text(hjust = 0.5))
plot(multi_logit1_plot)

## interaction term between Newspaper and Internet
multi_logit2 <- multinom(
  Knowledge_Class ~ FEMALE + GENERATION + 
    EDUCATION + INTEREST_RE + NEWSPAPER + ONLINE + NEWSPAPER * ONLINE, 
  data = voter2019_mlogit, reflevel = "Class1"
)
summary(multi_logit2)
multi_logit2_plot <- broom::tidy(multi_logit2, 
                                 conf.int = TRUE, 
                                 exponentiate = FALSE) %>% 
  dplyr::filter(term != "(Intercept)") %>% 
  mutate(term = recode(term, 
                       `FEMALE` = "Female", 
                       `GENERATION` = "Age",  
                       `EDUCATION` = "Education", 
                       `INTEREST_RE` = "Political Interest", 
                       `NEWSPAPER` = "Newspaper", 
                       `ONLINE` = "Internet", 
                       `NEWSPAPER:ONLINE` = "Newspaper * Internet")) %>% 
  mutate(term = factor(term, 
                       levels = rev(c("Female", "Age", "Education", 
                                      "Political Interest", "Newspaper", 
                                      "Internet", "Newspaper * Internet")))) %>% 
  ggplot(aes(x = term, y = estimate, 
             ymin = conf.low, ymax = conf.high)) + 
  geom_point() + geom_hline(yintercept = 0, color = "gray", linetype = "dashed") + 
  geom_pointrange() + coord_flip() + facet_wrap(~y.level, ncol = 3) + 
  labs(x = "Explanatory Variables", y = "Estimation", 
       title = "Multinomial Logistic Regression \nBaseline: Class1") + 
  theme(plot.title = element_text(hjust = 0.5))
plot(multi_logit2_plot)

## interaction term between Newspaper and TV
multi_logit3 <- multinom(
  Knowledge_Class ~ FEMALE + GENERATION + 
    EDUCATION + INTEREST_RE + TV + ONLINE + TV * ONLINE, 
  data = voter2019_mlogit, reflevel = "Class1"
)
summary(multi_logit3)
multi_logit3_plot <- broom::tidy(multi_logit3, 
                                 conf.int = TRUE, 
                                 exponentiate = FALSE) %>% 
  dplyr::filter(term != "(Intercept)") %>% 
  mutate(term = recode(term, 
                       `FEMALE` = "Female", 
                       `GENERATION` = "Age",  
                       `EDUCATION` = "Education", 
                       `INTEREST_RE` = "Political Interest", 
                       `TV` = "TV", 
                       `ONLINE` = "Internet", 
                       `TV:ONLINE` = "TV * Internet")) %>% 
  mutate(term = factor(term, 
                       levels = rev(c("Female", "Age", "Education", 
                                      "Political Interest", "TV", 
                                      "Internet", "TV * Internet")))) %>% 
  ggplot(aes(x = term, y = estimate, 
             ymin = conf.low, ymax = conf.high)) + 
  geom_point() + geom_hline(yintercept = 0, color = "gray", linetype = "dashed") + 
  geom_pointrange() + coord_flip() + facet_wrap(~y.level, ncol = 3) + 
  labs(x = "Explanatory Variables", y = "Estimation", 
       title = "Multinomial Logistic Regression \nBaseline: Class1") + 
  theme(plot.title = element_text(hjust = 0.5))
plot(multi_logit3_plot)
